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Abstract. The derivation of effective spin models describing the low 
energy magnetic properties of undoped CM02^planes is reinvestigated. Our 
study aims at a quantitative determination of the parameters of effective 
spin models from those of a multi-band model and is supposed to be rele- 
vant to the analysis of recent improved experimental data on the spin wave 
spectrum of La2Cu04. Starting from a conventional three-band model we 
determine the exchange couplings for the nearest and next-nearest neigh- 
bor Heisenberg exchange as well as for 4- and 6-spin exchange terms via a 
direct perturbation expansion up to 12th (14th for the 4-spin term) order 
with respect to the copper-oxygen hopping tp^. Our results demonstrate that 
this perturbation expansion does not converge for hopping parameters of the 
relevant size. Well behaved extrapolations of the couplings are derived, how- 
ever, in terms of Pade approximants. In order to check the significance of 
these results from the direct perturbation expansion we employ the Zhang- 
Rice reformulation of the three band model in terms of hybridizing oxygen 
Wannier orbitals centered at copper ion sites. In the Wannier notation the 
perturbation expansion is reorganized by an exact treatment of the strong 
site-diagonal hybridization. The perturbation expansion with respect to the 
weak intersite hybridizations is calculated up to 4th order for the Heisenberg 
coupling and up to 6th order for the 4-spin coupling. It shows excellent 
convergence and the results are in agreement with the Pade approximants 
of the direct expansion. The relevance of the 4-spin coupling as the leading 
correction to the nearest neighbor Heisenberg model is emphasized. 

Keywords: CM02-planes; three-band model; Heisenberg model; four spin 
interactions; quantum antiferromagnets. 



1 Introduction 



After the discovery of the high-Tc superconducting oxides it soon became 
clear that a minimum model for describing their electronic properties had to 
contain at least three bands derived from the copper crystal field state 
3dx2-y2 and from oxygen and 2py orbitals In a seminal paper Zhang 
and Rice showed |^ that the low energy physics of the three-band model is 
in fact contained in an effective single-band model, the type of model which 
was envisaged initially by Anderson 

The work presented here is concerned with a reinvestigation of the deriva- 
tion of effective single-band models from three-band models for Qt02-planes. 
The present paper will be confined to the study of undoped CM02-planes 
where the effective models contain spin degrees of freedom only. Effective 
low energy models are derived from high energy parent models via pertur- 
bative expansions 0. The focus of this work is placed on how to obtain 
high precision coupling constants for the effective spin models and what are 
the leading corrections to the familiar nearest neighbor Heisenberg model. 
For the system of strongly correlated electrons considered here the copper- 
oxygen hopping tpd is the expansion parameter of choice. The expansion in 
powers of tpd is, however, not straightforward due to its rather small radius 
of convergence. Therefore, expansions beyond the leading order are required 
for obtaining reliable results. This is probably the reason why in existing 
derivations of the magnetic Hamiltonian of CM02-planes the couplings are 
usually off by a factor of up to 2 0, . The dominant term in the effective 
Hamiltonian is the Heisenberg nearest neighbor exchange obtained in fourth 
order which is substantially corrected by higher order contributions which 
we will present up to twelfth order. In eighth order ring exchange processes 
start to contribute four-spin terms to the effective Hamiltonian which turn 
out to be not at all small 0. Our results are consistent with the recent 



interpretation [|T0[ of improved experimental data on the spin wave spectrum 
of La2Cu04 in terms of sizable four-spin exchange terms |TT|. In comparison 
to these four-spin terms second and third neighbor Heisenberg terms which 
also first appear in eighth order turn out to be rather tiny. We have cal- 
culated all these terms up to twelfth order (four-spin term up to fourteenth 
order) . It is evident from the results of these series expansions that physically 
relevant values of tpd are larger than the radius of convergence. We find, how- 
ever, that Fade approximants of these series expansions provide consistent 
extrapolations to the range of physically relevant model parameters. 
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There is an alternative approach to the perturbative treatment of three 
band models which shows a much better behavior of convergence and which 
we have also applied to obtain an independent check of the significance of 
the Pade approximants derived from the direct expansion. This approach has 
been introduced in the paper of Zhang and Rice on the three band model in 
which this model was reformulated in terms of hybridizing oxygen Wannier 
orbitals centered at the copper ion sites 0. In this notation the hopping 
Hamiltonian contains a large site-diagonal hybridization to which is easily 
treated exactly for each copper ion site and small intersite hybridizations 
which are then treated safely in a perturbative fashion. Along these lines 
Zhang and Rice achieved not only a clever rearrangement of the tpd pertur- 
bation series, but also a particularly transparent formulation of the physics of 
doped CM02-planes in terms of "spin" and "hole" states the latter of which 
are known as Zhang-Rice singlets. In the effective low energy model ("t-J 
model") obtained this way neighboring "spins" experience an exchange in- 
teraction J and "holes" interchange their position with neighboring "spins" 



via a hopping parameter t |T^. We will show that the leading contribution 



to the nearest neighbor Heisenberg exchange obtained in second order in the 
intersite hopping is sufficient to reproduce the major features found from the 
direct expansion up to realistic values of tpd, but is not sufficient for perfect 
agreement. Gided by sum rules for the hopping amplitudes in the Wannier 
representation we will then demonstrate how the agreement is systematically 
improved by including corrections of third and fourth order in the intersite 
hopping. The four-spin (up to sixth order) and further neighbor Heisenberg 
exchange terms will also be discussed in this context. 

The paper is organized as follows. In the following section the three band 
model used in this work is briefly reviewed together with its transformation 
into the Wannier representation. Section III describes the principles of the 
perturbative derivation of effective Hamiltonians as we will use it. Section 
IV is devoted to the direct expansion with respect to tpd and section V to 
the expansion in the Wannier representation. The results are summarized 
and conclusions are drawn in connection with the experimental evidence in 
section VI. 
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2 The three— band model 



In this section we will briefly present the three-band model from which our 
investigation is going to start and fix the notations used. For the purpose of 
this paper which is focusing on the feasibility of high precision determination 
of the parameters of effective spin models we will use a minimum three-band 
model with the Hamiltonian 

H = H, + Hu + (1) 

where the first term 

^ ^[^ddlad-l^a + ^pipl,\+n^/2,aPx,\+n^/2,a + Pl,\+ny/2,aPy,l+ny/2,cT)] (2) 

1,(T 

describes the energies of the 3d- and 2p-holes involved, the second term 

1 

describes the Coulomb repulsion of holes on the Cu^^ ions and the third 
term 

Hpd = tpdY.[4APx,l+n^/2,a+Py,l+ny/2,a-Px,l-n^/2,a~Py,l-ny/2,a)+^-C-]- (4) 
1,0- 

describes the hopping of holes between 3d- and neighboring 2p-sites. Copper 
3(i2:2_y2-orbitals are placed on a square lattice in the [x, |/)-plane which is 
spanned by unit vectors n^. and iiy and the vertices of which are labeled by 
the integer vector 1. Oxygen 2px- and 2py-orbitals are placed at the center 
of X- and y-bonds, respectively, between neighboring lattice sites. 

Typical parameters for the three-band model (|ip being used to model 
Cu02-planes are [13, H 

Apd = ep-ed = 3.6eV, U = 8eV, tpd = 1.3eV. (5) 

For a direct expansion with respect to the hopping parameter tp^ the Hamil- 
tonian (P is decomposed into 

H = Hl + V^ with Hl = H, + Hu and = Hp^. (6) 

Although the hopping amplitude tpd is smaller than the charge transfer energy 
Apd and than the Coulomb energy U it turns out that a direct expansion of 
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the parameters of an effective low energy model with respect to tpd, i.e. 
an expansion in powers of V^, does not work for the parameter set (^. 
We will demonstrate this later explicitly and we will estimate the radius of 
convergence of such a direct expansion as tp^ ~ U/16 = 0.5 eV. We will 
therefore work out this expansion to higher orders and will extract useful 
information from this expansion via Pade approximants. 

Zhang and Rice found an elegant way to reorganize the perturbation 
expansion by reformulating the three-band model in terms of hybridizing 
oxygen Wannier orbitals centered at the copper ion sites. The reformulated 
model is obtained after transforming the hopping term into momentum space 
representation using the Fourier transformed operators 

dl = 4^ E (7) 

and 

ri,+n./2,. = 4f E e-^'^('+-/^Vi,k,. (« = ^,y), (8) 
where L denotes the number of unit cells. With the form factor 

-f/i \ ■ ■ 2 , ■ 2 o /i COskx + COS ky , , 

/(k) = 2^sm' — + sm' y = 2^ 1 ^ (9) 

and the normalized hybridizing Wannier orbital in momentum space repre- 
sentation 

k k 
Wk,a = 2^(sin y ■ k,<x + sin -| ■ Py,k,^)//(k) (10) 

the hopping term reads 

Hp, = tp, E /(k)(4,.^k,. + <.4,.)- (11) 

k,cr 

Applying the Fourier transform (|^ to the Wannier operators wl_ ^ mutually 

orthogonal real space Wannier orbitals wl^ centered at the copper sites are 
obtained. In terms of these the hopping Hamiltonian finally takes the form 

Hpd = tpd E [Ti-mdlw^^^ + h.c], (12) 

l,m,(T 

where the Fourier coefficients 

TK^E/W^-^tl^/We- (13) 
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of the form factor have the full symmetry of the square lattice. Numerical 
values of these coefficients are given in Table 1. 



R 




(0,0) 


1.916183 


(±1,0),(0,±1) 


-0.280186 


(±1,±1) 


-0.047013 


(±2,0),(0,±2) 


-0.027450 


(±2,±1),(±1,±2) 


-0.013703 



Table 1. Numerical values for Tr 
The coefficients Tr satisfy the sum rules 



(/'(k)e^'^')k 

4 (1=(0,0)) 
-1 (l=(±l,0)or (0,±1)) 
(else). 



(14) 



which we are going to use later. Obviously, the site-diagonal amplitude ^(0,0) 
is much larger than all the other amplitudes and satisfies by itself the sum 
rule S(o,o) = 4 already to 91.8%. The amplitudes to the 4 first neighbors are 
almost 7 times smaller than T(o,o) and including them the sum rule S(o,o) = 4 is 
missed by only 0.35%. The amplitudes to further neighbors are much smaller 
again. One can show that in the limit of large distances the amplitudes drop 
asymptotically like 



Ti 



R 



{R 



00 



(15) 



27ri?3 

To write the Hamiltonian (0) also in terms of Wannier states non-hybridizing 
2p-orbitals orthogonal to the Wannier orbitals w have to be introduced. In 
momentum space representation they are given by 



2i(sin ^ ■ k,, 



k 



(16) 



and since the 2p-basis sets {px,Py) and {w,v) are unitarily equivalent one 
obtains 

He = T.[^ddlA,a + ^pHa^la + ^1%^)]- (17) 
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The Wannier representation in (|T^ and (|T^ allows a decomposition of the 
total Hamiltonian (|I]) into 

H = H^ + V'" (18) 

where a major part of the hopping term (^) is incorporated in the unper- 
turbed Hamiltonian. Using the shorthand notation 

^0 = 7'(o,o)^pd ~ 1-916 (19) 

the unperturbed Hamiltonian is chosen as Q 

1 

hi = Y^l^ddlJi^^ + epwl^w^^^ + toidl^w^^^ + wlji^^)] (20) 

a 

+f/(ii j^fij j^. 

The non-hybridizing orbital v can be ignored altogether in the minimum 
three-band model considered here since it is always completely filled. The 
local Hamiltonians hi act independently at each site 1. They are easily di- 
agonalized exactly. The perturbative part of the total Hamiltonian is then 
given by the intersite hopping terms in Wannier representation 



(21) 



l,m,< 



which are so small that they can safely be treated perturbatively for model 
parameters as given by (||). The separation of energy scales between and 
V"^ achieved through the use of the Wannier representation is so substantial 
that it is hard to understand why a controversy about the scenario proposed 
by Zhang and Rice 0] arose early on ||15|, |16|, |17|, [18], which was still quoted 
as unsettled in the review by Dagotto Leading order perturbative cal- 
culations using the Wannier representation were performed by many authors 
(see, e.g., M 



21, 22, 23, 24, 25, m 



3 Perturbative derivation of effective Hamil- 
tonians 

The perturbative derivation of effective Hamiltonians for correlated electron 
systems has a long history the early stages of which were summarized by 
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Takahashi in 1977 [^. In this paper Takahashi presents a particularly trans- 
parent description of the method and gives an explicit solution for the ef- 
fective Hamiltonian to arbitrary order. We will briefly recall Takahashi's 
approach here, because we are going to perform the perturbation expansions 
in this paper using his formulation and because we wish to avoid controversies 



about the proper use of the method like in [2i, 



It is assumed that the total Hamiltonian of a system is decomposed into 

H = Ho + V. (22) 

In the case of interest Hq has a degenerate subspace Uq of ground states with 
energy Eq. On switching on the perturbation V continuously the subspace 
Uq evolves continuously into the subspace U of the corresponding low energy 
eigenspace of H. Takahashi presents an explicit perturbative formula to all 
orders in V for an isometric linear transformation F: Uq ^ U describing the 
mapping of Uq onto U. In terms of F the effective Hamiltonian is then given 
by 

H,s = r^HT. (23) 

It acts in the subspace Uq of unperturbed eigenstates of Hq and has the 
same spectrum as the perturbed Hamiltonian H. In view of the explicit 
perturbation series of F it is a pure problem of book-keeping to set up the 
perturbation series for ifefr to any required order. In terms of the projection 
operator Pq onto the ground state subspace Uq and the resolvent operator 

S = (24) 

the full perturbation expansion up to fourth order is given by 

H,s = EoPo + PoVPo + PoVSVPo 

+P0VSVSVP0 - ^PoVPoVS^VPo - ^PoVS^VPoVPo 

+P0VSVSVSVP0 - ^PoVS^VPoVSVPo - ^PqVSVPoVS'^VPo 

+^PoVPoVPoVS^VPo + ^PqVS^VPoVPoVPo (25) 

-^PoVPoVS^VSVPo - ^PoVSVS^VPoVPo 

-^PoVPoVSVS^VPo - ^PoVS^VSVPoVPo. 
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For the purposes of the calculations in this paper we had to list this expansion 
up to twelfth order. The number of terms in the series grows exponentially 
with the order. To twelfth order the perturbation series contains 363721 
terms. 

For useful applications of the formal series of i^efr the unperturbed Hamil- 
tonian Hq has to be easily diagonalized such that matrix elements of the 
resolvent (p^ can be calculated explicitly. In this paper we will apply the 
perturbation expansion to the two Hamiltonian decompositions and ([18D 
were this condition on Hq is satisfied. We also will confine the analysis to 
undoped systems which implies that all terms in ifes containing PqVPq don't 
contribute. This reduces the number of twelfth order terms in ifcfr to 12341. 
In the direct expansion based on (|^) ground states can only be connected 
by an even number of hopping processes such that all terms with any odd 
number of V between two Pq don't contribute. This reduces the number of 
twelfth order terms to 3180. For the Wannier decomposition ([T8|) our analy- 
sis will be confined to sixth order. In this case, of the terms given in (pSf) 
only the second order term, the first third order term and the first three 
fourth order terms will contribute and up to sixth order 30 terms have to 
be taken into account. Notice that in the Wannier decomposition (|1^) odd 
order terms do contribute since Hq mixes d- and w-orbitals. 

4 Direct perturbation expansion 

In this section we are going to discuss the direct expansion with respect 
to tpd on the basis of the decomposition (^). In the undoped case that we 
are considering here the subspace of ground states Uq of the unperturbed 
Hamiltonian Hq contains all states without any p-holes and with a single 
(i-hole on each copper site. The effective Hamiltonian acting on Uq is thus 
a pure spin Hamiltonian acting on the spins 5 = 1/2 of the ci-sites. Due to 
the symmetry properties of the three band model (0) this Hamiltonian has 
got to be invariant under global spin rotations and under the space group 
of the square lattice. Only terms with an even number of spins are possible 
due to time reversal symmetry. The excited states of Hq are very simple and 
the excitation energies contain a Coulomb energy U for each rf-site with two 
holes and a charge transfer energy Ap^ for each p-hole. For a contribution of 
order n to the effective Hamiltonian one has to consider all sets of n hopping 
processes each of which defines a certain cluster of sites involved. Due to 
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the linked cluster theorem (which is bound to hold to keep the effective 
Hamiltonian extensive) only connected clusters are known to contribute. It 
is therefore sufficient to evaluate the various orders of Hcs on certain finite 
clusters. We have implemented the purely symbolic evaluation of the series 
expansion with a C-H- program. 

For simplicity we will disregard any constant energy shift in H^s since we 
want to focus on the effective spin Hamiltonian. The leading term in ifgff is 
then a fourth order nearest neighbor Heisenberg exchange Ji Si ■ S2 with the 



well known exchange coupling (see e.g. 



.(4) 4 2 44AU + \d) 

'''''' - AJ^2A;^ + - UAl, ■ ^26) 

In order to determine this coupling it is sufficient to calculate the amplitude 
of a spin flip process on a three-site cluster consisting of two neighboring 
d-sites and the p-site in between. In view of the identity Si ■ S2 = 5*1 + 
I (5*1" 5^ + S1S2) the coupling is given by twice the spin flip amplitude. 

In sixth order processes the six additional p-sites adjacent to the two 
d-sites can be visited by a hole. Therefore a nine site cluster would be 

(Pi) 

sufficient to calculate Jijjj.- Since in each individual exchange process at 
most one of the additional p-sites is visited the actual calculation can be 
confined to clusters of up to no more than four sites. Each of the six four- 
site clusters gives the same contribution to the sixth order coupling. In one 
such contribution either all four sites or only the three sites of the fourth 
order cluster will be involved in the exchange process. Therefore, the sixth 
order spin flip amplitude is given by six times the spin flip amplitude of 
the four-site cluster minus 5 times the spin flip amplitude of the three-site 
cluster. This type of reasoning would be dispensable in the sixth order case 
for which it was examplified here, but it is absolutely essential to make the 
higher order calculations feasible. It allows to reduce the maximum cluster 
size for the calculation of the nearest neighbor exchange from 17 to 8 in 
eighth order, from 31 to 9 in tenth order and from 43 to 12 in twelfth order. 

In eighth order ring exchange processes on an eight-site plaquette visiting 
four (i-sites are possible. These processes produce four-spin exchange terms 
in Hes- In cases where multi-spin terms are present the fewer-spin exchange 
terms can be inferred in the following way. Partial traces (i.e. traces over 
some of the spins) of any multi-spin term vanish. By forming the trace over 
some of the spins belonging to a cluster all exchange terms containing these 
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spins are therefore projected out. Applying this reasoning to the eight-site 
plaquette one obtains the two-spin exchange of a pair of spins by averaging 
over all configurations of the other spins contained in the plaquette. From 
time reversal invariance and hermiticity of Hes one can infer that the am- 
plitude of a spin flip process remains unchanged if all unfiipped spins of a 
cluster are inverted. This allows to reduce by a factor of 2 the number of 
configurations needed for the averaging. 

Along the lines described above we have calculated the nearest neighbor 
exchange coupling jjj. up to twelfth order in tpd- The full formula of the 
twelfth order result is given by Eq. ( [A.l| ) in Appendix A. Fig. 1 shows how 
the ratio ^jj./ Ji%,. varies with increasing tpd if the sixth, eighth, tenth 
and twelfth order terms are included (see the thick lines in Fig. 1). It is 
obvious from this plot that for the physically relevant values of tp^ as given 



1.2 
1 

0.8 

^1,dir 

|(4) 0.6 

^l.dir 

0.4 
0.2 




U=8 eV, Ap^=3.6 eV 



/ 

/' 




0.2 0.4 0.6 0.8 



6 th order . 
8 th order 
1 th order 
1 2 th order 
1-1 Pade 

1- 1 DIogP 

2- 2 Pade 



1.2 1.4 



Figure 1: Variation of ^i.dirZ-^i^ir with tp^ 



in Ji is smaller than simple estimates from would suggest, but 

it is also obvious that the radius of convergence of the direct perturbation 
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series considered here is much smaller than tp^ = 1.3 eV. The direct series 
determines Ji accurately only up to tpd ~ 0.5 eV. Extrapolations beyond the 
radius of convergence can be obtained, however, via Pade approximants of 
the series for ^^^/ We denote as m-n Pade the approximant with an 
j^th qj-^Iqj- numerator polynomial and an n**^ order denominator polynomial in 
the variable x = t^^. We have also constructed extrapolations via analogous 

Pade approximants for the logarithmic derivative of ^j^./ ^ii. which we 
denote by m-n DlogPade. The 1-1 and 2-2 Pades and the 1-1 DlogPade 
shown by the thin lines in Fig. 1 demonstrate the excellent convergence of 
this extrapolation procedure. The 0-1, 1-3 and 1-2 Pades and the 1-2 and 
2-1 DlogPades are not shown because they all differ from the 2-2 Pade by 
less than 3% for tpd < 1.3 eV and less than 4% for tp^ < 1.5 eV. From this 
observation we derive the estimate that they determine the nearest neighbor 
exchange coupling with an accuracy of better than 4%. Note the substantial 
reduction of the coupling in the range of physical interest, Ji = 0.33 dji. for 
tpd = 1.3 eV, in comparison to the lowest order result. 

The four-spin exchange terms which first appear in eighth order can be 
inferred from considering processes in which all four spins are flipped. Let 
us label the spins on the four d-sites of a square plaquette in cyclic order by 
numbers 1 to 4. There are 3 independent four-spin invariants, (Si-S2)(S3-S4), 
(S2 ■ S3)(S4 ■ Si) and (Si ■ S3)(S2 ■ S4) from which the four-spin exchange 
terms have to be formed. Due to the square point symmetry of our model 
the first two invariants always get the same exchange coupling in the effective 
Hamiltonian. This common coupling is given by twice the amplitude of the 
process which flips all spins of the initial state |1 t,2 |,3 1,4 |), since the 
third invariant doesn't contribute to this process. The exchange coupling of 
the third invariant can be inferred from considering an alternative four-spin 
flip process starting from the initial state 111,21,31,41). The sum of the 
couplings of the first invariant and the third invariant is given by four times 
the amplitude of this process. It turns out that this amplitude vanishes in 
eighth and tenth order. This implies that up to tenth order the four-spin 
exchange term has the form 

Jn [(Si ■ S2)(S3 ■ S4) + (S2 ■ S3)(S4 ■ Si) - (Si ■ S3)(S2 ■ S4)] (27) 

in analogy to what is known for the one band Hubbard model in fourth order 
i- 

The vanishing of the |1|,2|,3J.,4|) spin flip process up to tenth order can 
be easily understood as resulting from the linked cluster theorem, because for 
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these processes the plaquette (1,2,3,4) decomposes into two unhnked clusters, 
one of them containing the d-sites 1 and 2, the other containing sites 3 and 4. 
In twelfth order there are processes linking these two clusters and producing 
another four-spin term 

Jx (Si -83) (82-84) (28) 

in addition to (0). It has to be noted that in twelfth order clusters con- 
taining 6 (i-sites are created which produce six-spin terms in the effective 
Hamiltonian. In the calculation of the four-spin terms these six-spin terms 
have to be properly eliminated by the averaging procedure described above. 
The eighth order coupling constant of the four-spin term (^) is found to 

be 

,(8) _ 80t«,([/ + A,,) {U' + UA,, + Al,) 

Corrections up to fourteenth order are shown by Eq. ( |A.2|) in Appendix 
A together with the leading order contribution for Jx in Eq. ( |A.3| ). The 
variation of jj^/ ^ir with increasing tpd is shown in Fig. 2. Here, the 0- 
3 Fade (not shown) and the 1-2 and the 2-1 Fades as well as the 1-1 DlogFade 
seem to provide a rather accurate estimate with an uncertainty of about ±6% 
for tpd = 1.3 eV and an uncertainty of about ±15% for tpd = 1.5 eV. For 
tpd = 1.3 eV the coupling Ju is about 10 times smaller than suggested by the 
leading order term. We will consider the 1-1 DlogFade the most probable 

estimate of dir/^n dir* 

The leading contributions to second neighbor Heisenberg exchange terms 
like J2 8(0,0) " ^(1,1) to third neighbor Heisenberg exchange terms like 
J3 S(o,o) ■ S(o,2) are also obtained in eighth order. These couplings are given 
by 

Ul (uU^ + AU'A,d + 2U Al, + At 



and 



.(8) ^Pd^ ^pd^ ^pd, 



>3 



^(8) 4t^, {3U^ + 2U'Apd + 2UAl, + A^, 

f/3 A 



^3,dir - rr3 A7 ' y^^) 

^pd 

respectively. Corrections to these leading order expressions which we have 
calculated to twelfth order are given by Eqs. ( |A.4| ) and ( |A.5| ) in Appendix 
A. 

Figs. 3 and 4 show the variation of ^2dir/'^2dir '^dir/^sdir' respec- 
tively, with increasing tpd- The radii of convergence appear to be even smaller 
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U=8 eV, A ^=3.6 eV / 

pd I 




1 th order 
1 2 th order 
14 th order 

0- 2 Pade 

1- 2 Pade 

2- 1 Pade 
1-1 DIogP 



0.2 0.4 0.6 0.8 



Figure 2: Variation of Jn,dir/ -^cfdir with t 



pd- 



than in the case of Ji. The 1-1 Pade in Fig. 3 might indicate that J2 changes 
sign shghtly below tp^ = 1 eV, but the scattering of the various approximants 
doesn't allow definite conclusions on a change of sign. Since the 0-1 Dlog- 
Pade (not shown in Fig. 3) coincides to high precision with the 0-2 Pade we 
will consider this approximant as the most probable estimate for J2. The 
0-2 Pade for J3 shown in Fig. 4 turns upwards and has a pole at tpd ~ 2.1eV. 
The other three approximants shown appear to behave consistently and we 
will consider the 0-1 DlogPade as the most probable estimate for J3. Alto- 
gether, the Pade approximants for J2 and J3 scatter much more than those 
for Ji and Jn and provide less accurate estimates for J2 and J3. We do, 
however, learn from these extrapolations that for tpd = 1.3 eV both J2 and 
J3 also are reduced substantially in comparison to the leading order results 
( PPP and (PTD, J2 probably by a factor of as much as 10 and J3 probably by 
a factor of 5. As we will see later J2 and J3 are so small in absolute size that 
their accurate determination is less urgent for practical purposes. 
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Figure 3: Variation of J2,dir/<^2'dir with tp^. 

To describe the six-spin term resulting from twelfth order ring exchange 
processes on a double plaquette we label the six spins involved cyclically by 
numbers 1 to 6. Since the oxygen ion at the center of the double plaquette 
is not visited in twelfth order the six-spin term has the full symmetry of the 
hexagon formed by the six spins. The 15 independent invariants obtained by 
all pairings of the six spins into three scalar products group into the 5 
operators with hexagonal symmetry Oi to O5 given by Eq. ( |A.6| ) in Appendix 
A. With the same type of arguments which led to ( p7|) we conclude that the 
twelfth (and fourteenth) order six-spin term has the form 

^□(01 + 02-03 + 04-05). (32) 
r(i2) 



The exchange coupling given by Eq. ( |A.7| ) in Appendix A was calculated 



from ring exchange processes which flip all spins of the state |1|,2|,3|,4J, 
,5T,6i). 

For a comparison of the relative sizes of the various couplings we first show 
in Fig. 5 the leading perturbative contributions of all couplings determined, 
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Figure 4: Variation of Jz,dir/ J3 dir with tp^. 

in units of In the range of physically relevant model parameters the 

four-spin coupling Jn is by far the largest correction to the nearest neighbor 
two-spin coupling Ji . The second and third neighbor Heisenberg couplings J2 
and J3 are much smaller and are in fact comparable to the six-spin coupling 
Jq. This scenario agrees with what is known from perturbation expansions 
for the single band Hubbard model ^ and from cluster calculations for 
the three band model 0. 

A quantitative comparison of the best approximants for the various cou- 
plings with Ji (represented by its 2-2 Fade) is shown in Fig. 6 where we 
have denoted the m-n Fade for the coupling Jj by Jj [m, ri]. For the model 
parameters Jn is almost one order of magnitude smaller than Ji and 
the couplings J2 and J3 are almost another order of magnitude smaller. The 
four-spin coupling Jn therefore has to be considered an important modifi- 
cation of the simple nearest neighbor Heisenberg model, whereas the second 
and third neighbor Heisenberg couplings J2 and J3 (as well as the six-spin 
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Figure 5: Comparison of the leading order terms of the various couphngs. 



couphng J^) may be ignored as correction at the level of about 3%. 



5 Expansion in the Wannier representation 

In this section we are going to discuss the alternative perturbation expan- 
sion based on the decomposition (|TB|) of the three band Hamiltonian. The 
unperturbed Hamiltonian pO|) consists of independent local Hamiltonians hi 
for each site which are easily diagonalized. In the one hole sector the local 
Hamiltonian has two S = 1/2 eigenstates where mutually orthogonal linear 
combinations of d-hole and w-hole orbitals are occupied. The two hole sec- 
tor contains a non-hybridized 5* = 1 triplet and three hybridized singlets the 
lowest one of which is the Zhang-Rice singlet. In the three hole sector again 
two S = 1/2 doublets are found. The four hole sector and the zero hole sector 
each contain one trivial S = state. The lower doublet in the one hole sector 
acts as the local ground state doublet of the undoped system. All the other 
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Figure 6: Comparison of the best Fade approximants. 



states will show up as intermediate excited states at sufficiently high orders 
of the expansion with respect to the perturbation (pTj). We have diagonal- 
ized the local Hamiltonian numerically. A simple analytic formula cannot 
be obtained in the general case since for the three singlets in the two hole 
sector a (3x3)-matrix has to be diagonalized. Simple analytic expressions 
for the solution in this sector would be available only in the symmetric case 
^pd = U/2. The perturbation expansion with respect to (plf ) was performed 
using a combination of symbolic and numerical routines. 

It is instructive to analyse the radius of convergence t^^ of the local Hamil- 
tonian hi which does depend on tpd via (19). This radius of convergence can 
be determined from studying the branch points of the eigenvalues of the lo- 
cal Hamiltonian in the complex tp^-plane. Without going into any details 
we wish to summarize this analysis here by stating that tp^ = 0.469 eV for 
the values of U and Apd given in (|). This value agrees well with the val- 
ues estimated from the Fade approximants. This is not surprising since the 
expansion with respect to the small perturbation (|21| ) is expected to con- 
verge well and should not much modify t^^ as defined above from the local 
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Hamiltonian. 
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Figure 7: Variation of Ji,w/Ji,kr with tpd- 

In what follows we will plot the variation of the various coupling constants 
with the hopping tpd in analogy to the presentations in the previous figures 
by measuring all couplings in units of their lowest order term in the direct 
expansion (if not otherwise stated). Fig. 7 shows our results for the nearest 
neighbor exchange Ji^w In the present context the leading contribution to 
Ji^yj is obtained from the simple second order hopping process described by 
the term PqVSVPq of (pS]). This second order contribution is depicted by the 
thick dotted line in Fig. 7. It is satisfying that this simple second order result 
reproduces quite nicely the decrease of ^i/^^tlir with increasing tpd as given 
by the Pade approximants of Fig. 1. On the other hand, there is, however, 
a systematic deviation in the overall size of the coupling; even for small tpd 
the coupling Jf ,^ is too large by about 15%. The discrepancy at small tpd 
is largely reduced by taking into account the third order terms derived from 
PqVSVSVPq and, finally, jf^ is in satisfying agreement with the 2-2 Pade 
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Figure 8: Variation of Ju^JJ^ir with t,,. 
of the direct expansion. The deviation of Ji^2 from for small tpd is 

(2) 

explained quantitatively by re-expanding J| to second order with respect 
to to- Referring to (0) we obtain 

ji% ~ (2r(o,o)7'(i,o))^^i% {tpd -* 0) (33) 

which explains the 15% deviation because (2T(o,o)7'(i,o))^ = 1-153. Since J^^^^ 
collects all fourth order terms it has to coincide with after re-expanding 
it to tp^. How this happens becomes particularly clear if one looks at the 
sum rule S(i^o) of (0)- This sum rule states that 2T(o,o)7'(i,o) + ^(i,o) = ~1 if 
we denote by r(i^o) the sum of all terms in S(i,o) (infinitely many) which don't 
contain T(o,o)- Squaring this sum rule we obtain the relation 

(2T(o,o)T(i,o))' + 2(2T(o,o)T(i,o))r(i,o) + r^^o) = 1 (34) 

from which we can read off the contributions of various orders of the Wannier 
expansion to J[%t. in the limit of small tp^. The first term represents the 
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contribution of J},^ discussed above. The second term contains only one 
factor of ^(0,0) and results from third order terms in which due to r(i q) = 
0.073775 exhaust the relation ( |34D to 1 — r^^ = 0.994557; this explains 

why Ji^l is slightly smaller than in the limit of small tpd (see Fig. 7). 
Finally, the term r^^ comes from fourth order terms in V"^ which contribute 
only about 0.5% for small tpd but change sign and get more important as tpd 
increases. 

In our calculation of third and fourth order contributions to Ji^^ we have 
made extensive use of the sum rule S(i^o)- In third order terms the exchange 
path for a spin flip process involves an arbitrary third copper ion site whose 
spin is not flipped. The sum of the spin flip amplitudes over all these third 
sites contains a lattice sum which is simply r(i o). The calculation of the 
fourth order is more involved since one has to discriminate between spin flip 
processes which don't visit another site and those which visit one or two 
more sites. The lattice sums appearing in this order cannot be completely 
determined from sum rules, but sum rules considerably simplify their calcu- 
lation. Four-spin terms which here appear in fourth order are eliminated by 
averaging as described in the previous section. 

Results for the four-spin coupling Jn^w are shown in Fig. 8. The leading 
fourth order contribution again nicely reproduces qualitatively the decrease 
with increasing tpd known from Fig. 2. After the above discussion the devia- 
tion observed in the small tpd limit is not surprising. In fact, a quantitative 
understanding of this deviation follows from looking at the fourth power of 
the S(i,o) sum rule: (2T(o,o)T(i,o) +'^(1,0))'' = 1- With (2T(o,o)T(i,o))^ = 1.329 we 
understand why jj^]^ is about 33% too large for small tpd- The substantial re- 
duction of the deviation by the fifth order contributions are also understood 
quantitatively from the identity (2T(o,o)r(i^o))'' + 4(2T(o,o)7'(i,o))'^'^(i,o) = 0.964 
(see Fig. 8). Including the sixth order terms we find the almost negligible de- 
viation of (2T(o,o)T(i,o))" + 4(2r(o,o)T(i,o))V(i,o) + 6(2T(o,o)T(i,o))V2^ o) = 1.0017 
and the overall agreement with the 1-1 DlogPade is quite satisfying. 

The analysis of the second and third neighbor exchange J2 and J3 is 
more complicated in the Wannier representation since there are low order 
contributions which have to be cancelled completely by higher order terms 
before results of any significance emerge. We therefore show these couplings 
in Figs. 9 and 10 not in units of their eighth order counterparts from the 
direct expansion but in units of It is quite obvious that for any lattice 

vector 1 there is a second order contribution J^^ to the Heisenberg coupling 
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Figure 9: Variation of J2,«)/<^itiir with tp^ 



between two spins separated by 1 which in analogy to (|3^) behaves hke 



(2) 



(2T(o,o)Ti)^jj^jij. 



{tpd 0). 



(35) 



The cancellation of this contribution by higher order terms is understood by 
invoking the sum rule 2T(o,o)ri + ^"1 = for further neighbors which squared 
gives the relation 



(2T(o,o)Ti)2 + 2(2T(o,o)Ti)ri + rf = 0. 



The numbers (2T(o,o)7'( 



0.0325 and (2r(o,o)T(2,o))' 



(36) 



0.0111 coin- 



(2) (2) 

cide perfectly with the behavior of J2^^ and J^,^ for small tpd as shown in 
Figs. 9 and 10. We also understand from (^) why the inclusion of the third 
order doesn't reduce the deviation from but just changes its sign. In fourth 
order Fig. 9 shows that in agreement with (p6D the terms proportional to t^^ 
in J2^w vanish. This is, however, only a partial solution of the cancellation 
problem since there are still terms proportional to t^^ which according to 
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Figure 10: Variation of Js,^/ J^^jjj. with tpd- 

Fig. 9 even have the wrong sign and cancellation of which would only be 
achieved by extending the expansion to sixth order. For tpd > 0.8 eV 
the third and fourth order results shown in Figs. 9 and 10 at least have the 
right sign and the same order of magnitude as the Fade estimates from the 
previous section. We have to conclude that the accurate determination of the 
further neighbor couplings J2 and J3 in the Wannier representation would be 
very demanding. This points at definite limitations of this approach. 

6 Conclusions 



In the present paper we have discussed the derivation of high precision effec- 
tive spin Hamiltonians for the low energy sector of a three band model for 
Crt02~pla'nes. By two methods we have demonstrated that it is possible to 
overcome the convergence problems of the tpd perturbation series. Using the 
direct expansion with respect to tpd we have derived precise values for the 
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most important couplings via Pade approximants. The direct expansion has 
the advantage of a particularly simple unperturbed Hamiltonian and a very 
nicely localized perturbative Hamiltonian which makes high order symbolic 
expansions feasible. Using the Wannier representation we have confirmed 
the results from the direct expansion by a method with much better conver- 
gence properties. The expansion in the Wannier representation is, however, 
rendered more difficult by a more complicated unperturbed Hamiltonian and 
a less well localized perturbative Hamiltonian and by the necessity of a non- 
symbolic (i.e. numerical) series expansion. We have also shown that for 
precise values of the coupling constants the leading orders of the Wannier 
expansion are not sufficient. 

The work in the present paper was confined to the most rudimentary 
three band model since our main goal was to demonstrate the feasibility of 
the derivation of accurate effective Hamiltonians. Nevertheless we will use 
our results here for a fit of the couplings Ji = 151.9 meV and Jn/ Ji = 0.24 ex- 
tracted recently from a fit to the experimental dispersion of La2CuOi |11] us- 
ing self consistent spin-wave theory [jTO[. Assuming the typical (though some- 
what arbitrary) model parameters U = 8eV and Ap^ = 3.6 eV from (|^) we 
obtain from the value Ji = 151.9 meV the estimate 1.422 eV < tpd < 1.454 eV 
for the hopping parameter of the minimum model showing an uncertainty in 

of 2% due to the uncertainty of our Pade extrapolations. With tpd in this 
range our estimate for Jn results in 0.19 < Ju/Ji < 0.25 which is in good 
agreement with the result from fTO . 



For proper applications to cuprate materials this work will have to be 
extended to more realistic three band models including, in particular, a direct 



oxygen-oxygen hopping tpp UTS] , . The relevance of four-spin exchange has 



been stressed also for the two-leg ladder system LaQCagCu240 41 to which 
the analysis presented here can be applied as well. 
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Appendix A 



This Appendix contains the more voluminous formulae from the direct ex- 
pansion of section IV. These formulae can be easily used to derive the Pade 
approximants discussed in section IV. 

With (^6|) the Taylor series for the nearest neighbor exchange coupling is 



dir 



t{4) 



I -it 



4 (5 ?7 + 2 A 



Vd) 



801f/3 + 164 t/^A 



+ 



pd 



12 A3, 



8505 + 9602 Apd + 908 A^, - 240 U A^, 



48 AJ, 



+ 



U'AIAU + A,,f 
t^rf(758199 + 1587453 Ap^ + 890808 Aj, + 52603 f/^ A^, 
-6611 U' AJ, + 4566 U' A^ + 2559 U A^, + 483 A^ / 
(4f/^Aj,(f/ + A,,f) + 0(Ol. (A.l) 



The series for the four-spin coupling with the leading contribution 
given by 



IS 



□ ,dir 



j: 



(8) 



□ ,dir 



ti 



4 (llf/3 + 14f/2Apd + 8f/A2, + 2A 



^pd 



+ 



Al, {U + Ap,) {u^ + UAp, + A 
ijrf(56569 [/^ + 161892 Apd + 168480 f/^ Aj, + 76092 [/^ A: 
+9096 A J, - 7008 Aj, - 3960 f/ Ajj, - 792 A^J/ 
(40 Al, {U + A,,)3 + f/ A,, + AJJ) 
-tL(410565 f/^ + 1487797 Apd + 2034672 A 



3 



pd 



+1264452 A^, + 296152 Aj, - 48240 U'^ A J, 

-49264 f/^Aj, - 13464 f/Aj, - 1584 A^/ 

(10 AJ, {U + A,,)4 (f/2 + f/ A,, + AJ,)) + 0(tj,) 

The leading contribution to the four-spin coupling 



(A.2) 



IS 



J. 



(12) 



2 ( (489 f/^ + 1016 Apd - 72 Aj, - 2232 f/^ A 



pd 

3392 f/3 Aj, - 2784 Aj, - 1280 U A^ - 256 Ajj)/ 



(f/^AS(t/ + A 



(A.3) 
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With ( PUD and (^) the series for the second and third neighbor two-spin 
couphngs are given by 



'^2Air ~ '^2,dir 



4 (142 + 169 A,, + 36 Ag, + 10 [/ A^, + 2 A^,) 
A^,([/ + A,,)(llf/3 + 4f/2A,, + 2f/A^, + A3,) 

tjrf (82083 f/^ + 171784 Apd + 99154 Aj^ + 10848 Aj^ 
+1420 f/=^ Aj, + 290 U' Al, + 120 U Aj, + 24 Ajj/ 
(4 f/2 AJ, (f/ + A,,r (11 1/^ + 4 A,, + 2 f/ AJ, + Aj,)) 
+0(tyl (A.4) 



and 



7 - 

•^3,dir ~ "^S.dir 



1 



2(72^^ + 94^3A,, + 39[/^A^, + 20[/A^, + 4A^J 
AlAU + A,,)i3U^ + 2U^ A,, + 2U Al, + A^,) 
tj^(47947 f/^ + 111156 Apd + 90704 f/^ Aj^ + 43130 f/^ A; 
+21424 Ajrf + 8174 Aj^ + 2632 U A^^ + 440 Ajj/ 
(8 Al, (U + Apdf (3 + 2 f/2 Ap, + 2 f/ Aj, + AjJ) 
+0(ty]. (A.5) 

The 5 six-spin invariants for a hexagonal plaquette are 

01 = (Si-S2)(S3-S4)(S5-S6) + (S2-S3)(S4-S5)(S6-Si) 

02 = (Si-S4)(S2-S6)(S3-S5) + (S2-S5)(S3-Si)(S4-S6) 

+ (S3-S6)(S4-S2)(S5-Si) 

03 = (Si-S4)(S2-S5)(S3-S6) 

04 = (Si-S2)(S3-S6)(S4-S5) + (S2-S3)(S4-Si)(S5-S6) 

+ (S3-S4)(S5-S2)(S6-Si) 

05 = (Si-S2)(S3-S5)(S4-S6) + (S2-S3)(S4-S6)(S5-Si) 

+ (83 ■ S4)(S5 ■ Si)(S6 ■ S2) + (S4 ■ S5)(S6 ■ S2)(Si • S3) 

+(85 ■ S6)(Si ■ S3)(S2 ■ S4) + (Se ■ Si)(S2 ■ S4)(S3 ■ S5). (A.6) 
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The leading contribution to the six-spin couphng in Eq. (B^) is found to be 



(12) ^ 336 Q {U + A,,) (3 [/^ + 6 A,, + 8 A^, + 6 ^ Ag, + 3 A^,) 

(A.7) 
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